Cost of facilitation: invasive grasses limit recruitment of benefactor shrubs.

ecoblender

alex.filazzola

Abstract

Shrubs facilitate the abundance and productive of annual plants in desert ecosystems.However, these shrub microhabitats favour plant species with competitive life histories. Ephedra californica is a dominant shrub in the San Joaquin desert that has been identified as a facilitator. Here, we explore the factors that limit the Ephedra californica recruitment into the San Joaquin desert including substrate, water availability, and herbivory. We also explore the role of the invasive grass Bromus madritensis on limiting establishment of E. californica. Vegetation surveys were conducted in the field during 2013 and collected seed was then used to conduct two greenhouse trials. The first explore germination and establishment techniques for E. california and the optimal substrate. The second examined E. californica establishment in present of the invasive B. madritensis responding to different water levels and herbivory. These results can have implications for land managers in the San Joaquin Valley to maintain native shrub biodiversity in the region.

---

## load libraries
library(tidyverse)
library(OIsurv)
library(lsmeans)

## load data
substrate <-read.csv("data/Ephedra.substrate.csv")
recruit <-read.csv("data/Ephedra.recruitment.csv")
landscape <-read.csv("data/ephedra.landscape.csv")

##load functions
## inverse hyperbolic sine transformation 
ihs <- function(x) {
    y <- log(x + sqrt(x ^ 2 + 1))
    return(y)
}
se <- function(x, ...) sqrt(var(na.omit(x))/length(na.omit(x)))

Ephedra recruitment at landscape

## compare ephedra size and density across landscape
landscape[,"area.group"] <- as.numeric(cut(landscape$Log.area,20))
hist(landscape$area.group)

mean.area <- landscape %>% group_by(area.group) %>%  summarize(rdm=mean(RDM.2013), area=mean(Log.area), density=mean(Shrub.density))

ggplot(mean.area) + geom_point(aes(x=density, y=area), size=3)+ylab("log (area of shrub)")+ theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),axis.line = element_line(colour = "black"),panel.background = element_blank(),axis.text=element_text(size=14),axis.title=element_text(size=16)) + xlab(expression("shrub density (m"^"2"*")")) +  stat_smooth(method="lm", formula= y~x,aes(x=density, y=area))

m1 <- lm(area~density, data=mean.area)
summary(m1)

## 
## Call:
## lm(formula = area ~ density, data = mean.area)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.82510 -0.26068  0.01197  0.22216  1.16732 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)     9.3959     0.9852   9.537 1.84e-08 ***
## density     -1706.0323   223.6888  -7.627 4.81e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.5569 on 18 degrees of freedom
## Multiple R-squared:  0.7637, Adjusted R-squared:  0.7506 
## F-statistic: 58.17 on 1 and 18 DF,  p-value: 4.806e-07

ggplot(mean.area) + geom_point(aes(x=density, y=rdm), size=mean.area$area+1)+ylab("residual dry matter")+ theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),axis.line = element_line(colour = "black"),panel.background = element_blank(),axis.text=element_text(size=14),axis.title=element_text(size=16)) + xlab(expression("shrub density (m"^"2"*")")) +  stat_smooth(method="lm", formula= y~x,aes(x=density, y=rdm))

m2 <- lm(rdm~density, data=mean.area)
summary(m2)

## 
## Call:
## lm(formula = rdm ~ density, data = mean.area)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.9038 -0.8124 -0.3260  0.2697  4.4984 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    16.431      2.634   6.239 6.94e-06 ***
## density     -2081.306    597.979  -3.481  0.00267 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.489 on 18 degrees of freedom
## Multiple R-squared:  0.4023, Adjusted R-squared:  0.3691 
## F-statistic: 12.11 on 1 and 18 DF,  p-value: 0.002669

Ephedra optimal substrate

substrate[,"Sand"] <- as.factor(substrate$Sand)
substrate[,"census"] <- as.factor(substrate$census)

m1 <- aov(height ~ Micro * Sand,  data=subset(substrate, census==10))
summary(m2) ## nothing significant

## 
## Call:
## lm(formula = rdm ~ density, data = mean.area)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.9038 -0.8124 -0.3260  0.2697  4.4984 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    16.431      2.634   6.239 6.94e-06 ***
## density     -2081.306    597.979  -3.481  0.00267 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.489 on 18 degrees of freedom
## Multiple R-squared:  0.4023, Adjusted R-squared:  0.3691 
## F-statistic: 12.11 on 1 and 18 DF,  p-value: 0.002669

## compare shade on survival of ephedra
my.surv <- Surv(as.numeric(substrate$census), substrate$survival)
fit1 <- survfit(my.surv~Micro, data=substrate)
summary(fit1)

## Call: survfit(formula = my.surv ~ Micro, data = substrate)
## 
##                 Micro=Open 
##  time n.risk n.event survival std.err lower 95% CI upper 95% CI
##     1    960      11    0.989 0.00343        0.982        0.995
##     2    880      29    0.956 0.00681        0.943        0.969
##     3    800      41    0.907 0.00987        0.888        0.927
##     4    720      42    0.854 0.01221        0.830        0.878
##     5    640      67    0.765 0.01504        0.736        0.795
##     6    481      31    0.715 0.01647        0.684        0.748
##     7    400      62    0.604 0.01901        0.568        0.643
##     8    241      31    0.527 0.02108        0.487        0.570
##     9    160      31    0.425 0.02366        0.381        0.474
##    10     80      31    0.260 0.02730        0.212        0.320
## 
##                 Micro=Shade 
##  time n.risk n.event survival std.err lower 95% CI upper 95% CI
##     1    960      33    0.966 0.00588        0.954        0.977
##     2    880      43    0.918 0.00897        0.901        0.936
##     3    800      42    0.870 0.01117        0.849        0.892
##     4    720      40    0.822 0.01290        0.797        0.848
##     5    640      34    0.778 0.01422        0.751        0.807
##     6    480      17    0.751 0.01521        0.721        0.781
##     7    400      31    0.692 0.01725        0.659        0.727
##     8    240      16    0.646 0.01958        0.609        0.686
##     9    160      15    0.586 0.02317        0.542        0.633
##    10     80      13    0.491 0.03099        0.433        0.555

par(mar=c(4.5,4.5,.5,.5))
plot(fit1, col="white", ylim=c(0.2,1), xlim=c(0.9,10.1), ylab="Estimated Survival Function", xlab="Census", cex.lab=1.5, cex.axis=1.3)
lines(1:10,summary(fit1)[[10]][1:10], col="#E69F00", lty=2) ## upper
lines(1:10,summary(fit1)[[6]][1:10], col="#E69F00", lwd=2) ## value sun
lines(1:10,summary(fit1)[[11]][1:10], col="#E69F00", lty=2) ## lower
lines(1:10,summary(fit1)[[10]][11:20], col="#56B4E9", lty=2) ## upper
lines(1:10,summary(fit1)[[6]][11:20], col="#56B4E9", lwd=2) ## value shade
lines(1:10,summary(fit1)[[11]][11:20], col="#56B4E9", lty=2) ## lower
legend(1.8, 0.35, c("Sun","Shade"), lty=1, lwd=3, col=c("#E69F00","#56B4E9"), cex=1.5)

coxph.fit <- coxph(my.surv ~ Micro, method="breslow", data=substrate)
coxph.fit 

## Call:
## coxph(formula = my.surv ~ Micro, data = substrate, method = "breslow")
## 
##               coef exp(coef) se(coef)     z       p
## MicroShade -0.2802    0.7557   0.0786 -3.56 0.00037
## 
## Likelihood ratio test=12.8  on 1 df, p=0.000342
## n= 1920, number of events= 660

## compare sand on survival of ephedra
my.surv <- Surv(as.numeric(substrate$census), substrate$survival)
fit2 <- survfit(my.surv~Sand, data=substrate)
summary(fit2)

## Call: survfit(formula = my.surv ~ Sand, data = substrate)
## 
##                 Sand=0 
##  time n.risk n.event survival std.err lower 95% CI upper 95% CI
##     1    384       7    0.982 0.00683        0.968        0.995
##     2    352      11    0.951 0.01125        0.929        0.973
##     3    320      13    0.912 0.01506        0.883        0.942
##     4    288      12    0.874 0.01799        0.840        0.910
##     5    256      14    0.827 0.02106        0.786        0.869
##     6    192       6    0.801 0.02289        0.757        0.847
##     7    160      12    0.741 0.02695        0.690        0.795
##     8     96       6    0.694 0.03120        0.636        0.758
##     9     64       6    0.629 0.03794        0.559        0.708
##    10     32       6    0.511 0.05325        0.417        0.627
## 
##                 Sand=25 
##  time n.risk n.event survival std.err lower 95% CI upper 95% CI
##     1    384      10    0.974 0.00813        0.958        0.990
##     2    352      16    0.930 0.01331        0.904        0.956
##     3    320      20    0.872 0.01772        0.838        0.907
##     4    288      18    0.817 0.02075        0.777        0.859
##     5    256      15    0.769 0.02292        0.726        0.815
##     6    192       7    0.741 0.02441        0.695        0.791
##     7    160      14    0.676 0.02776        0.624        0.733
##     8     96       7    0.627 0.03137        0.568        0.692
##     9     64       7    0.558 0.03714        0.490        0.636
##    10     32       7    0.436 0.05007        0.348        0.546
## 
##                 Sand=50 
##  time n.risk n.event survival std.err lower 95% CI upper 95% CI
##     1    384      11    0.971 0.00851        0.955        0.988
##     2    352      16    0.927 0.01350        0.901        0.954
##     3    320      19    0.872 0.01765        0.838        0.907
##     4    288      20    0.812 0.02098        0.771        0.854
##     5    256      18    0.755 0.02343        0.710        0.802
##     6    192       9    0.719 0.02512        0.672        0.770
##     7    160      17    0.643 0.02848        0.589        0.701
##     8     96       8    0.589 0.03178        0.530        0.655
##     9     64       8    0.516 0.03697        0.448        0.593
##    10     32       8    0.387 0.04823        0.303        0.494
## 
##                 Sand=75 
##  time n.risk n.event survival std.err lower 95% CI upper 95% CI
##     1    384       9    0.977 0.00772        0.962        0.992
##     2    352      18    0.927 0.01361        0.900        0.954
##     3    320      16    0.880 0.01716        0.847        0.915
##     4    288      17    0.828 0.02025        0.790        0.869
##     5    256      28    0.738 0.02422        0.692        0.787
##     6    192      13    0.688 0.02624        0.638        0.741
##     7    160      26    0.576 0.02976        0.521        0.637
##     8     96      13    0.498 0.03266        0.438        0.566
##     9     64      13    0.397 0.03612        0.332        0.474
##    10     32      13    0.236 0.04058        0.168        0.330
## 
##                 Sand=100 
##  time n.risk n.event survival std.err lower 95% CI upper 95% CI
##     1    384       7    0.982 0.00683        0.968        0.995
##     2    352      11    0.951 0.01125        0.929        0.973
##     3    320      15    0.907 0.01554        0.877        0.937
##     4    288      15    0.859 0.01891        0.823        0.897
##     5    256      26    0.772 0.02349        0.727        0.819
##     6    193      13    0.720 0.02596        0.671        0.773
##     7    160      24    0.612 0.03000        0.556        0.674
##     8     97      13    0.530 0.03351        0.468        0.600
##     9     64      12    0.431 0.03755        0.363        0.511
##    10     32      10    0.296 0.04372        0.222        0.395

test <- substrate %>% group_by(census,Sand) %>%  summarize(count=sum(survival),avg=mean(survival))
test <- data.frame(test)

par(mar=c(4.5,4.5,.5,.5))
plot(fit2, col="white", ylim=c(0.2,1), xlim=c(0.9,10.1), ylab="Estimated Survival Function", xlab="Census", cex.lab=1.5, cex.axis=1.3)
lines(1:10,summary(fit2)[[10]][1:10], col="#E69F00", lty=2) ## upper
lines(1:10,summary(fit2)[[6]][1:10], col="#E69F00", lwd=2) ## value sand 0
lines(1:10,summary(fit2)[[11]][1:10], col="#E69F00", lty=2) ## lower
lines(1:10,summary(fit2)[[10]][11:20], col="#56B4E9", lty=2) ## upper
lines(1:10,summary(fit2)[[6]][11:20], col="#56B4E9", lwd=2) ## value sand 25
lines(1:10,summary(fit2)[[11]][11:20], col="#56B4E9", lty=2) ## lower
lines(1:10,summary(fit2)[[10]][21:30], col="#009E73", lty=2) ## upper
lines(1:10,summary(fit2)[[6]][21:30], col="#009E73", lwd=2) ## value sand 50
lines(1:10,summary(fit2)[[11]][21:30], col="#009E73", lty=2) ## lower
lines(1:10,summary(fit2)[[10]][31:40], col="#CC79A7", lty=2) ## upper
lines(1:10,summary(fit2)[[6]][31:40], col="#CC79A7", lwd=2) ## value sand 75
lines(1:10,summary(fit2)[[11]][31:40], col="#CC79A7", lty=2) ## lower
lines(1:10,summary(fit2)[[10]][41:50], col="#D55E00", lty=2) ## upper
lines(1:10,summary(fit2)[[6]][41:50], col="#D55E00", lwd=2) ## value sand 100
lines(1:10,summary(fit2)[[11]][41:50], col="#D55E00", lty=2) ## lower
legend(1.7, 0.44, c("0%","25%","50%","75%","100%"), lty=1, lwd=3, col=c("#E69F00","#56B4E9","#009E73","#CC79A7","#D55E00"), cex=1.2, title="Sand")

coxph.fit <- coxph(my.surv ~ Sand, method="breslow", data=substrate)
coxph.fit 

## Call:
## coxph(formula = my.surv ~ Sand, data = substrate, method = "breslow")
## 
##          coef exp(coef) se(coef)    z       p
## Sand25  0.263     1.301    0.138 1.91  0.0563
## Sand50  0.365     1.441    0.135 2.71  0.0068
## Sand75  0.579     1.785    0.130 4.47 7.7e-06
## Sand100 0.450     1.568    0.133 3.39  0.0007
## 
## Likelihood ratio test=23.3  on 4 df, p=0.000109
## n= 1920, number of events= 660

## raw number of germinants end trials

end <- subset(substrate, census==10)

## compare microsite

end.means <- end %>% group_by(Micro) %>% summarize(eph=mean(survival), eph.se=se(survival), Height=mean(height, na.rm=T), Height.se=se(height)) %>% gather(measure, value, 2:5) %>% separate(measure, c("Eph",".se")) ## calculate metrics and put long format

## Warning: package 'bindrcpp' was built under R version 3.4.2

## Warning: Too few values at 4 locations: 1, 2, 5, 6

end.means[is.na(end.means$.se),3] <- "mean"
end.means <- end.means %>% spread(.se, value)
end.means <- data.frame(end.means)

end.means[,"lower"] <- end.means[,"mean"]-end.means[,"se"]
end.means[,"upper"] <- end.means[,"mean"]+end.means[,"se"]

survival <- subset(end.means, Eph=="eph")

plot1 <- ggplot(survival, aes(x=Micro, y=mean)) + geom_bar(stat="identity") +
  geom_errorbar(aes(ymin=lower, ymax=upper), width=.2)+ylab("average surviving Ephedra plant")+ theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),axis.line = element_line(colour = "black"),panel.background = element_blank(),axis.text=element_text(size=14),axis.title=element_text(size=16)) + xlab("Microsite") 

height <- subset(end.means, Eph=="Height")

plot2 <- ggplot(height, aes(x=Micro, y=mean)) + geom_bar(stat="identity") +
  geom_errorbar(aes(ymin=lower, ymax=upper), width=.2)+ylab("average height of Ephedra plant")+ theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),axis.line = element_line(colour = "black"),panel.background = element_blank(),axis.text=element_text(size=14),axis.title=element_text(size=16)) + xlab("Microsite") 

require(gridExtra)

## Loading required package: gridExtra

## 
## Attaching package: 'gridExtra'

## The following object is masked from 'package:dplyr':
## 
##     combine

grid.arrange(plot1, plot2, ncol=2)

m1 <- glm(survival~ Sand* Micro, data=end, family="binomial")
anova(m1, test="Chisq")

## Analysis of Deviance Table
## 
## Model: binomial, link: logit
## 
## Response: survival
## 
## Terms added sequentially (first to last)
## 
## 
##            Df Deviance Resid. Df Resid. Dev Pr(>Chi)   
## NULL                         159     188.21            
## Sand        4   4.7398       155     183.47 0.315056   
## Micro       1  10.7288       154     172.75 0.001055 **
## Sand:Micro  4   8.6588       150     164.09 0.070217 . 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

lsmeans(m1, pairwise~Sand)

## NOTE: Results may be misleading due to involvement in interactions

## $lsmeans
##  Sand     lsmean          SE df   asymp.LCL    asymp.UCL
##  0    -1.4663371   0.4529108 NA   -2.354026  -0.57864819
##  25   -1.3671838   0.4643107 NA   -2.277216  -0.45715153
##  50   -8.7830342 494.5226042 NA -978.029528 960.46345952
##  75   -0.4236489   0.3831780 NA   -1.174664   0.32736619
##  100  -0.8047190   0.3872983 NA   -1.563810  -0.04562817
## 
## Results are averaged over the levels of: Micro 
## Results are given on the logit (not the response) scale. 
## Confidence level used: 0.95 
## 
## $contrasts
##  contrast    estimate          SE df z.ratio p.value
##  0 - 25   -0.09915331   0.6486236 NA  -0.153  0.9999
##  0 - 50    7.31669717 494.5228116 NA   0.015  1.0000
##  0 - 75   -1.04268814   0.5932568 NA  -1.758  0.3987
##  0 - 100  -0.66161811   0.5959263 NA  -1.110  0.8012
##  25 - 50   7.41585049 494.5228222 NA   0.015  1.0000
##  25 - 75  -0.94353482   0.6020048 NA  -1.567  0.5185
##  25 - 100 -0.56246480   0.6046358 NA  -0.930  0.8852
##  50 - 75  -8.35938531 494.5227526 NA  -0.017  1.0000
##  50 - 100 -7.97831529 494.5227559 NA  -0.016  1.0000
##  75 - 100  0.38107003   0.5448168 NA   0.699  0.9567
## 
## Results are averaged over the levels of: Micro 
## Results are given on the log odds ratio (not the response) scale. 
## P value adjustment: tukey method for comparing a family of 5 estimates

## survival
end.means <- end %>% group_by(Sand) %>% summarize(eph=mean(survival), eph.se=se(survival))
end.means <- data.frame(end.means)

end.means[,"lower"] <- end.means[,"eph"]-end.means[,"eph.se"]
end.means[,"upper"] <- end.means[,"eph"]+end.means[,"eph.se"]

ggplot(end.means, aes(x=Sand, y=eph)) + geom_bar(stat="identity", fill="#56B4E9") +
  geom_errorbar(aes(ymin=lower, ymax=upper), width=.2)+ylab("average surviving Ephedra plant")+ theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),axis.line = element_line(colour = "black"),panel.background = element_blank(),axis.text=element_text(size=14),axis.title=element_text(size=16)) + xlab("Percentage of sand in soil") 

## height

end.means <- end %>% group_by(Sand) %>% summarize(eph=mean(height, na.rm=T), eph.se=se(height, na.rm=T))
end.means <- data.frame(end.means)

end.means[,"lower"] <- end.means[,"eph"]-end.means[,"eph.se"]
end.means[,"upper"] <- end.means[,"eph"]+end.means[,"eph.se"]

ggplot(end.means, aes(x=Sand, y=eph)) + geom_bar(stat="identity", fill="#56B4E9") +
  geom_errorbar(aes(ymin=lower, ymax=upper), width=.2)+ylab("average surviving Ephedra plant")+ theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),axis.line = element_line(colour = "black"),panel.background = element_blank(),axis.text=element_text(size=14),axis.title=element_text(size=16)) + xlab("Percentage of sand in soil") 

## most germinants
maximum.germ <- substrate %>% group_by(Micro, Sand, Rep) %>% summarize(eph.max=max(survival))

max.eph <- maximum.germ %>% group_by(Sand) %>% summarize(eph=mean(eph.max,na.rm=T),eph.se=se(eph.max, na.rm=T)) 
max.eph <- data.frame(max.eph)

max.eph[,"lower"] <- max.eph[,"eph"]-max.eph[,"eph.se"]
max.eph[,"upper"] <- max.eph[,"eph"]+max.eph[,"eph.se"]


ggplot(max.eph, aes(x=Sand, y=eph)) + geom_point(fill="black", size=4) +
  geom_errorbar(aes(ymin=lower, ymax=upper), width=.2)+ylab("total Ephedra germinants")+ theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),axis.line = element_line(colour = "black"),panel.background = element_blank(),axis.text=element_text(size=14),axis.title=element_text(size=16)) + xlab("Percentage of sand in soil") + ylim(0,1) + stat_smooth(method="lm", formula= y~poly(x,2),data=max.eph, aes(x=as.numeric(Sand), y=eph), lwd=1, color="black", lty=2)

m1 <- lm(eph~poly(as.numeric(Sand),2), data=max.eph)
summary(m1)

## 
## Call:
## lm(formula = eph ~ poly(as.numeric(Sand), 2), data = max.eph)
## 
## Residuals:
##        1        2        3        4        5 
## -0.01250  0.01875  0.01875 -0.04375  0.01875 
## 
## Coefficients:
##                            Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                 0.60625    0.01768  34.295 0.000849 ***
## poly(as.numeric(Sand), 2)1  0.09882    0.03953   2.500 0.129612    
## poly(as.numeric(Sand), 2)2 -0.23385    0.03953  -5.916 0.027402 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.03953 on 2 degrees of freedom
## Multiple R-squared:  0.9538, Adjusted R-squared:  0.9075 
## F-statistic: 20.62 on 2 and 2 DF,  p-value: 0.04624

Test of limiting factors - brome

recruit[is.na(recruit)] <- 0

## water
water <- subset(recruit, Treatment == "control" | Treatment=="water")

## germination
water.germ <- water %>% group_by(Density, Lvl) %>% summarize(eph=mean(ephedra.end ),brome=mean(brome.begin))

ggplot(water.germ) + geom_jitter(aes(x=Density, y=brome, fill=Lvl, color=Lvl), size=2)+ylab("number of Bromus germinants")+ theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),axis.line = element_line(colour = "black"),panel.background = element_blank(),axis.text=element_text(size=14),axis.title=element_text(size=16))

m1 <- glm(brome.begin~Density *Lvl, data=water, family=poisson)
anova(m1, test="Chisq")

## Analysis of Deviance Table
## 
## Model: poisson, link: log
## 
## Response: brome.begin
## 
## Terms added sequentially (first to last)
## 
## 
##             Df Deviance Resid. Df Resid. Dev Pr(>Chi)    
## NULL                          299    1275.05             
## Density      1   849.52       298     425.53   <2e-16 ***
## Lvl          2     4.50       296     421.03   0.1053    
## Density:Lvl  2     0.00       294     421.03   0.9985    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## biomass
water.bio <- water %>% group_by(Density, Lvl) %>% summarize(eph=mean(Ephedra.biomass ),brome=mean(brome.biomass))

ggplot(water.bio) + geom_jitter(aes(x=Density, y=brome, fill=Lvl, color=Lvl), size=2)+ylab("final biomass of Brome")+ theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),axis.line = element_line(colour = "black"),panel.background = element_blank(),axis.text=element_text(size=14),axis.title=element_text(size=16))

## clipping
clip <- subset(recruit, Treatment == "control" | Treatment=="clipped")

## germination
clip.germ <- clip %>% group_by(Density, Lvl) %>% summarize(eph=mean(ephedra.end ),brome=mean(brome.begin))

ggplot(clip.germ) + geom_jitter(aes(x=Density, y=brome, fill=Lvl, color=Lvl), size=2)+ylab("number of Bromus germinants")+ theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),axis.line = element_line(colour = "black"),panel.background = element_blank(),axis.text=element_text(size=14),axis.title=element_text(size=16))

## biomass
clip.bio <- clip %>% group_by(Density, Lvl) %>% summarize(eph=mean(Ephedra.biomass ),brome=mean(brome.biomass))

ggplot(clip.bio) + geom_jitter(aes(x=Density, y=brome, fill=Lvl, color=Lvl), size=2)+ylab("final biomass of Brome")+ theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),axis.line = element_line(colour = "black"),panel.background = element_blank(),axis.text=element_text(size=14),axis.title=element_text(size=16))

m1 <- glm(brome.begin~Density *Lvl, data=clip, family=poisson)
anova(m1, test="Chisq")

## Analysis of Deviance Table
## 
## Model: poisson, link: log
## 
## Response: brome.begin
## 
## Terms added sequentially (first to last)
## 
## 
##             Df Deviance Resid. Df Resid. Dev Pr(>Chi)    
## NULL                          299    1332.93             
## Density      1   863.30       298     469.63   <2e-16 ***
## Lvl          2     3.66       296     465.97   0.1605    
## Density:Lvl  2     0.04       294     465.93   0.9814    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## shade
shade <- subset(recruit, Treatment == "control" | Treatment=="shade")

## germination
shade.germ <- shade %>% group_by(Density, Lvl) %>% summarize(eph=mean(ephedra.end ),brome=mean(brome.begin))

ggplot(shade.germ) + geom_jitter(aes(x=Density, y=brome, fill=Lvl, color=Lvl), size=2)+ylab("number of Bromus germinants")+ theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),axis.line = element_line(colour = "black"),panel.background = element_blank(),axis.text=element_text(size=14),axis.title=element_text(size=16))

## biomass
shade.bio <- shade %>% group_by(Density, Lvl) %>% summarize(eph=mean(Ephedra.biomass ),brome=mean(brome.biomass))

ggplot(shade.bio) + geom_jitter(aes(x=Density, y=brome, fill=Lvl, color=Lvl), size=2)+ylab("final biomass of Brome")+ theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),axis.line = element_line(colour = "black"),panel.background = element_blank(),axis.text=element_text(size=14),axis.title=element_text(size=16))

m1 <- glm(brome.begin~Density *Lvl, data=shade, family=poisson)
anova(m1, test="Chisq")

## Analysis of Deviance Table
## 
## Model: poisson, link: log
## 
## Response: brome.begin
## 
## Terms added sequentially (first to last)
## 
## 
##             Df Deviance Resid. Df Resid. Dev  Pr(>Chi)    
## NULL                          299    1221.52              
## Density      1   804.97       298     416.55 < 2.2e-16 ***
## Lvl          2    13.34       296     403.21  0.001267 ** 
## Density:Lvl  2     0.22       294     402.98  0.894507    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Test of limiting factors - Ephedra

## water
## germination

ggplot(water.germ) + geom_jitter(aes(x=Density, y=eph, fill=Lvl, color=Lvl), size=2)+ylab("number of Ephedra germinants")+ theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),axis.line = element_line(colour = "black"),panel.background = element_blank(),axis.text=element_text(size=14),axis.title=element_text(size=16))

m1 <- glm(ephedra.end~Density *Lvl, data=water, family=poisson)
anova(m1, test="Chisq")

## Analysis of Deviance Table
## 
## Model: poisson, link: log
## 
## Response: ephedra.end
## 
## Terms added sequentially (first to last)
## 
## 
##             Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL                          299     220.66         
## Density      1  1.76764       298     218.90   0.1837
## Lvl          2  0.12445       296     218.77   0.9397
## Density:Lvl  2  0.05766       294     218.72   0.9716

## biomass

ggplot(water.bio) + geom_jitter(aes(x=Density, y=eph, fill=Lvl, color=Lvl), size=2)+ylab("final biomass of Ephedra")+ theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),axis.line = element_line(colour = "black"),panel.background = element_blank(),axis.text=element_text(size=14),axis.title=element_text(size=16))

## clipping

## germination

ggplot(clip.germ) + geom_jitter(aes(x=Density, y=eph, fill=Lvl, color=Lvl), size=2)+ylab("number of Ephedra germinants")+ theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),axis.line = element_line(colour = "black"),panel.background = element_blank(),axis.text=element_text(size=14),axis.title=element_text(size=16))

## biomass

ggplot(clip.bio) + geom_jitter(aes(x=Density, y=eph, fill=Lvl, color=Lvl), size=2)+ylab("final biomass of Ephedra")+ theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),axis.line = element_line(colour = "black"),panel.background = element_blank(),axis.text=element_text(size=14),axis.title=element_text(size=16))

m1 <- glm(ephedra.end~Density *Lvl, data=clip, family=poisson)
anova(m1, test="Chisq")

## Analysis of Deviance Table
## 
## Model: poisson, link: log
## 
## Response: ephedra.end
## 
## Terms added sequentially (first to last)
## 
## 
##             Df Deviance Resid. Df Resid. Dev Pr(>Chi)  
## NULL                          299     220.71           
## Density      1  2.73641       298     217.97  0.09808 .
## Lvl          2  0.34733       296     217.63  0.84058  
## Density:Lvl  2  0.47792       294     217.15  0.78745  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## shade

## germination

ggplot(shade.germ) + geom_jitter(aes(x=Density, y=eph, fill=Lvl, color=Lvl), size=2)+ylab("number of Ephedra germinants")+ theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),axis.line = element_line(colour = "black"),panel.background = element_blank(),axis.text=element_text(size=14),axis.title=element_text(size=16))

##biomass
ggplot(shade.bio) + geom_jitter(aes(x=Density, y=eph, fill=Lvl, color=Lvl), size=2)+ylab("final biomass of Ephedra")+ theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),axis.line = element_line(colour = "black"),panel.background = element_blank(),axis.text=element_text(size=14),axis.title=element_text(size=16))

m1 <- glm(brome.begin~Density *Lvl, data=shade, family=poisson)
anova(m1, test="Chisq")

## Analysis of Deviance Table
## 
## Model: poisson, link: log
## 
## Response: brome.begin
## 
## Terms added sequentially (first to last)
## 
## 
##             Df Deviance Resid. Df Resid. Dev  Pr(>Chi)    
## NULL                          299    1221.52              
## Density      1   804.97       298     416.55 < 2.2e-16 ***
## Lvl          2    13.34       296     403.21  0.001267 ** 
## Density:Lvl  2     0.22       294     402.98  0.894507    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1